#/bin/bash
usage="
Author  : Hyunmin Kim (hyun.kim@ucdenver.edu)
Tool    : calculate density from a bam file 
Usage   : $0 [options] <bam> 

[options]:
    -a <int> : shifting size (default 100)
    -n       : normalize by RPM (read per million)
"
SHIFT=100;
NORM=0;
while getopts "hna:" arg; do
    case $arg in
        a) SHIFT=$OPTARG;;
        n) NORM=1;;
        *) echo "$usage"; exit 0;;
    esac
done
shift $(( OPTIND - 1 ));
if [ $# -ne 1 ];then echo "$usage"; exit 0; fi
FBAM=$1
if [ ! -f $FBAM.idx ];then
    samtools index $FBAM
fi
TOT=0;
if [ $NORM == 1 ];then
    TOT=$( samtools idxstats $FBAM | awk 'BEGIN{ tot=0; }{ tot += $3; } END{ print tot;}' )
fi
## for each chromosome generate count files
GINFO=( `samtools idxstats $FBAM | awk '$1 != "*"'` )
for (( i=0; i< ${#GINFO[@]}; i+=4))
do
    CHROM=${GINFO[$i]}; LEN=${GINFO[$i+1]};
    echo "processing $CHROM .. " >&2
    samtools view -bq 1 $FBAM $CHROM | bamToBed \
    | awk -vSHIFT=$SHIFT  -vLEN=$LEN '{OFS="\t";
        chrom=$1;start=$2;end=$3;strand=$6; 
        if(strand == "+"){
            pos = start + int(SHIFT);
        }else{
            pos = end - int(SHIFT) - 1;
        }       
        if(pos < 0){ pos = 0;
        }else if(pos > LEN - 1){ pos = LEN-1; }
        print pos;
    }' \
    | sort -nk1 | uniq -c | awk -vCHROM=$CHROM -vTOT=$TOT 'BEGIN{ OFS="\t";FS=" ";};{ count=$1; if(TOT > 0){ count=count*1000000/TOT} print CHROM,$2,$2+1,count;}'
done
